

%==============================================================
% Driver program to calculate the V_coinc and V_coherent for a 
% network of gravitational wave detectors.
%==============================================================

clear all;

                   
  % make a source sky grid with good resolution:
  angularResolution = 0.05;
  skyPositions = sinusoidalMap(angularResolution);

  % extract the spherical angles associated with the sky grid:
  thetas = skyPositions(:,1);
  phis = skyPositions(:,2);

  detNames = {'LHO','LLO','Virgo','INDIGO'};
  nDetectors = length(detNames);

  psi = 0.0;
  iota = 0.0;
  mTot = 2.8; %double Neutron star binary
  indigoDetPsi = 0.0; %degrees

  [rmax_coher] = Rmax_coherent_comb(iota, psi, phis, thetas, detNames, indigoDetPsi, mTot);

  dOmega = skyPositions(:, 4);

  coher_vol.HLV = (1/3) * sum ((rmax_coher.HLV).^3 .* dOmega);

  coher_Reff.HLV = (3/(4*pi)) * (coher_vol.HLV)^(1/3);

  coher_vol.HLVI = (1/3) * sum ((rmax_coher.HLVI).^3 .* dOmega)

  coher_Reff.HLVI = (3/(4*pi)) * (coher_vol.HLVI)^(1/3)

  coher_f_ratio = (rmax_coher.HLVI - rmax_coher.HLV)./ rmax_coher.HLV;


  %========== Plotting code ============
  try 
      Plot_Rmax_Coherent;
  catch
      fprintf(2, 'Could not execute Plot_Rmax_Coherent - Could be License issue ?\n');
  end 


  %====================================================================================
  %====================================================================================
  k = 1;
  indigoDetPsi = 0; %degrees
  % Vary the detector's arm geometry and plot the change in Rmax
  while indigoDetPsi <= 90
       
       clear rmax_coher;  
       [rmax_coher] = Rmax_coherent_comb(iota, psi, phis, thetas, detNames, indigoDetPsi, mTot);

       mean_coher_rmax_HLVI(k) = mean(rmax_coher.HLVI);

       k = k+1;
       indigoDetPsi = indigoDetPsi + 5;
  end
  figure(367); clf;
  plot([0:5:90], mean_coher_rmax_HLVI, 'b-', 'LineWidth', 2);
  title('Avg. Coherent Distance vs. IndIGO Arm Geometry -- aaaa', 'FontSize', 12);
  xlabel('IndIGO Arm Geometry / Deg', 'FontSize', 14);
  ylabel('Avg. Coher. Dist. / Mpc', 'FontSize', 14);
  axis tight; box on;
  set(gca, 'FontSize', 12);
  print('-dpng', 'Coherent_VaryDetPsi_aaaa.png');

  %================================================================================================
  %================================================================================================
 
  k = 1;
  indigoDetPsi = 0; %degrees
  mTot = 2.0; %solar masses 
  clear mean_coher_rmax_HLVI;
  % Vary the masses and plot the change in Rmax
  while mTot <= 200
       
       clear rmax_coher;  
       [rmax_coher] = Rmax_coherent_comb(iota, psi, phis, thetas, detNames, indigoDetPsi, mTot);
        
       mean_coher_rmax_HLVI(k) = mean(rmax_coher.HLVI);
       mean_coher_rmax_HLV(k) = mean(rmax_coher.HLV);

       k = k+1;
       mTot = mTot + 10;
  end
  figure(467); clf;
  plot([2:10:200], mean_coher_rmax_HLVI, 'b-', 'LineWidth', 2); 
  hold on;
  plot([2:10:200], mean_coher_rmax_HLV, 'b-o', 'LineWidth', 2); 
  hold off;
  title('Average Coherent Distance vs. Total Mass of Binary System-- aaaa', 'FontSize', 12);
  xlabel('Mass_{Total}/ M_{sun}', 'FontSize', 14);
  ylabel('Avg. Coher. Dist. / Mpc', 'FontSize', 14);
  axis tight; box on;
  set(gca, 'FontSize', 12);
  legend('(HLVI)', '(HLV)');
  print('-dpng', 'Coherent_VaryMtot_aaaa.png');
 	

 
   
  





